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Abstract. We report new radial velocity observations of GP Vel / HD 77581, the optical companion to the 
eclipsing X-ray pulsar Vela X-l. Using data spanning more than two complete orbits of the system, we detect 
evidence for tidally induced non-radial oscillations on the surface of GP Vel, apparent as peaks in the power 
spectrum of the residuals to the radial velocity curve fit. By removing the effect of these oscillations (to first 
order) and binning the radial velocities, we have determined the semi-amplitude of the radial velocity curve of 
GP Vel to be K a — 22.6 ± 1.5 km s" 1 . Given the accurately measured semi-amplitude of the pulsar's orbit, 
the mass ratio of the system is 0.081 ± 0.005. We are able to set upper and lower limits on the masses of the 
component stars as follows. Assuming GP Vel fills its Roche lobe then the inclination angle of the system, i, is 
70.1° ± 2.6°. In this case we obtain the masses of the two stars as M x = 2.27 ± 0.17 Mq for the neutron star and 
M = 27.9 ± 1.3 Mq for GP Vel. Conversely, assuming the inclination angle is i = 90°, the ratio of the radius of 
GP Vel to the radius of its Roche lobe is /3 = 0.89±0.03 and the masses of the two stars are M x = 1.88±0.13 Mq 
and M — 23.1 ± 0.2 Mq. A range of solutions between these two sets of limits is also possible, corresponding to 
other combinations of i and 0. In addition, we note that if the zero phase of the radial velocity curve is allowed 
as a free parameter, rather than constrained by the X-ray ephemeris, a significantly improved fit is obtained with 
an amplitude of 21.2 ± 0.7 km s _1 and a phase shift of 0.033 ± 0.007 in true anomaly. The apparent shift in the 
zero phase of the radial velocity curve may indicate the presence of an additional radial velocity component at 
the orbital period. This may be another manifestation of the tidally induced non-radial oscillations and provides 
an additional source of uncertainty in the determination of the orbital radial velocity amplitude. 
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1. Introduction 

Binary systems containing an X-ray pulsar are very im- 
portant astrophysically as they can offer a direct measure- 
ment of the neutron star mass. The mass ratio of the sys- 
tem, q, is simply given by the ratio of the radial velocity 
semi-amplitudes for each component: 



M n 



Ka 



(1) 



where M x is the neutron star mass, M is the mass of 
the optical component, K Q is the semi- amplitude of the 
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radial velocity of the optical component, and K x is the 
semi-amplitude of the radial velocity of the neutron star. 
In addition, for an elliptical orbit, it can be shown that, 
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where i is the inclination of the orbital plane to the 
line of sight, e is the eccentricity and P is the period of 
the orbit. 
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We therefore have a means of calculating the mass of 
the neutron star if the orbits of the two components and 
the inclination of the system are known. Such a combi- 
nation is possible in an eclipsing X-ray binary system, in 
which the neutron star is a pulsar. X-ray pulse timing de- 
lays around the neutron star orbit yield the value of K x , 
and conventional radial velocity measurements from opti- 
cal spectra yield K . 

A value for i can be obtained from the following ap- 
proximations: 



Table 1. Measured amplitudes of the radial velocity curve 
of GP Vol 
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where Rl is the radius of the optical companion's 
Roche lobe, (3 = R /Rl is the ratio of the radius of the 
optical companion to that of its Roche lobe, o' is the sep- 
aration between the centres of masses of the two compo- 
nents, and 9 e is the eclipse half-angle. A, B, and C have 
been determined by Rappaport & Joss (1983) to be: 



A « 0.398 - 0.026ft 2 + 0.004ft 3 
B « -0.264 + 0.052ft 2 - 0.015ft 3 
C « -0.023 - 0.005ft 2 



(6) 
(7) 
(8) 



where ft is the ratio of the rotational period of the 
companion star to the orbital period of the system. Note 
that, whilst R^/a! is expected to be constant for a given 
system, R\ J1 a' and (3 will vary around the orbit for a 
system which has an appreciable eccentricity. 

Only seven eclipsing X-ray binary pulsars are known 
(namely Her X-l, Cen X-3, Vela X-l, SMC X-l, LMC X-4, 
QV Nor and OA01657-415). Orbital parameters for the 
first six of these are still relatively poorly known, whilst 
the counterpart to OA01657-415 has only recently been 
identified (Chakrabarty et al 2002) and no optical radial 
velocity curve has been measured. In addition, the mass 
of the neutron star in a seventh eclipsing X-ray binary 
(4U1700-37) has recently been determined by Clark et 
al (2002). This system does not contain an X-ray pulsar 
though, and the mass determination is based on a Monte- 
Carlo modelling method which relics on the spectral type 
of the companion and is thus highly uncertain. Vela X-l 
is the only one of these systems to have an eccentric orbit, 
and apart from OAO 1657-415 is the one with the longest 
orbital period. For these reasons, and others discussed be- 
low, an accurate determination of the mass of the neutron 
star in Vela X-l has always been difficult to obtain. 

2. The history of Vela X-l 

Vela X-l was first detected by a rocket-borne experiment 
(Chodil et al 1967) and subsequent observations (see, for 
example, Giacconi ct al 1972), suggested that the source 



K 1 km s- 1 


Reference 


37-45 


Wallerstein (1974) 


26 ±0.7 


Zuiderwijk 1974 


20 ±1 


van Paradijs et al (1976) 


21.75 ±1.15 


van Paradijs et al (1977) 


21.8 ± 1.2 


Rappaport & Joss (1983) 


18.0-28.2 


(95% conf. range) van Kerkwijk et al (1995) 


17.8 ±1.6 


Stickland et al (1997) 


22 


(Correction to the above) Barziv et al (2001) 


21.7 ±1.6 


Barziv et al (2001) 


22.6 ± 1.5 


this paper 



was highly variable. Using data from the OS 0-7 satellite, 
Ulmer et al (1972) demonstrated evidence for intensity 
variations which were interpreted as eclipses with a pe- 
riod of around 9 days. An optical counterpart, GP Vel / 
HD 77581 (a B0. 5 giant with m v = 6.8) was identified by 
Brucato & Kristian (1972) and Hiltner et al (1972), based 
on its ultra-violet excess and radial velocity variations. 

An X-ray pulse period of 283 s was subsequently 
discovered using the SAS-3 satellite (Rappaport & 
McClintock 1975, McClintock et al 1976). Timing obser- 
vations of these pulses by Rappaport, Joss & McClintock 
(1976) allowed the radial velocity semi-amplitude of the X- 
ray component to be measured as K x = 273 ± 9 km s _1 , 
and the eccentricity of the system, e <~ 0.1. Later work, 
using Hakucho and Tenma data (Deeter et al 1987), gave 
a value of a x sinz = 113.0 ± 0.4 light seconds for the pro- 
jected semi-major axis. Using values for the orbital period 
P and the eccentricity of the orbit e, the corresponding 
K x value may be calculated according to 



2tt a x sin i 

p"(l_ e 2)l/2 



(9) 



as K x — 276±1 km s _1 . The most recent, and accurate, 
values determined from X-ray pulse timing analysis are 
those obtained with BATSE on board CGRO reported 
by Bildsten et al (1997). They obtain a value a x sini = 
113.89 ± 0.13 light seconds which corresponds to K x = 
278.1 ± 0.3 km s _1 when combined with their accurate 
values for the orbital period, P = 8.964368 ± 0.000040 d, 
and eccentricity, e = 0.0898 ± 0.0012. 

The X-ray eclipse duration appears to be quite vari- 
able, and somewhat energy dependent. For example, 
Forman et al (1973) obtained a value for the half angle of 
the eclipse of 9 e = 38° ±1° using the Uhuru satellite and 
Charles et al (1976) obtained a value of 9 e = 39°. 8 ± 0°.4 
using the Copernicus satellite, whereas Watson & Griffiths 
(1977) quote a value of 9 e = 33°. 8 ± 1°.3 obtained using 
Ariel V data. However the earlier two experiments were 
at much softer energies than Watson & Griffiths observed, 
and softer X-rays are much more likely to be absorbed by 
circumstellar material, thus extending the observed eclipse 
time. 
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Early determinations of K Q were made by Wallerstein 
(1974) and Zuiderwijk et al (1974), see Table 1. As regu- 
lar X-ray pulsations had yet to be discovered at this point, 
assumptions had to be made about the mass of the opti- 
cal component in order to estimate the mass of the com- 
pact object. Zuiderwijk et al (1974) obtained a value of 
M x > 2.5 ± 0.3 M Q , and suggested that such a large mass 
coupled with the lack of regular pulsations indicated that 
the compact object was a black hole. The discovery of reg- 
ular X-ray pulsations provided a means of determining the 
masses of both components directly, and also ruled out the 
possibility of the compact object being a black hole. 

Van Paradijs et al (1976) combined their K Q value ob- 
tained from 26 coude spectrograms (Table 1) with the 
K x value of Rappaport & McClintock (1975). Using X- 
ray eclipse data, they determined that i > 74°, and thus 
arrived at a mass of 1.6 ± 0.3 M Q for the neutron star. 
Van Paradijs et al (1977) subsequently refined their K Q 
value using yet more photographic spectra (Table 1) and 
Rappaport & Joss (1983) revised K further (Table 1) 
obtaining a neutron star mass estimate of 1.85lo1jo M© 
by combining data from a number of sources, includ- 
ing Watson & Griffiths (1977), and Rappaport, Joss & 
Stothcrs (1980), and performing a Monte Carlo analysis 
to estimate the uncertainties. 

More recently, van Kerkwijk et al (1995) made further 
optical observations of GP Vel, and discovered strong de- 
viations from a pure Keplerian velocity curve, which were 
auto-correlated within a single night, but not from one 
night to another. It was suggested that the variable grav- 
itational force exerted by the neutron star as it travels 
around its eccentric orbit excites short-lived oscillations 
on the surface of the optical component which affect the 
measured radial velocity. From their K value (Table 1) 
van Kerkwijk et al (1995) obtained M x = 1.9±°l M Q . A 
significantly lower value for K Q (Table 1) was obtained 
from observations using the IUE satellite by Stickland, 
Lloyd & Radziun-Woodham (1997). However, Barziv et 
al (2001) report that the analysis of these IUE data was 
subject to an error and a correct analysis yields a value 
consistent with those previously measured (Table 1) thus 
solving the discrepancy. 

The most recent measurement of the optical radial 
velocity curve of GP Vel (Barziv et al 2001) made use 
of 183 spectra obtained over a nine month campaign in 
order to try to average out the deviations reported by 
van Kerkwijk et al (1995). Although they were quite suc- 
cessful in averaging out these excursions, they were left 
with different, phase- locked deviations in the radial veloc- 
ity curve. Despite this they determined an accurate K 
value (Table 1) and set a limit on the neutron star mass 
of M x sin 3 i = 1.78 ± 0.15 M . 

3. Observations 

To take account of the deviations from Keplerian motion 
that are present in the radial velocity curve of GP Vel, we 
adopted a somewhat different approach to that employed 



by Barziv et al (2001). Instead of trying to average the 
velocity excursions over many orbital periods, we chose 
to carry out a comprehensive radial velocity study with 
maximum phase coverage over several consecutive orbital 
cycles. In this way we aimed to track the velocity excur- 
sions closely and thus model and ultimately remove them. 

The observations of GP Vel reported here were 
made using the 74-inch telescope at Mount Stromlo 
Observatory, Canberra, Australia. In order to cover two 
complete orbits, twenty-one continuous nights were ob- 
tained between 1996 February 1 and 1996 February 21. 
Only two nights were completely lost (the 9th and 17th) 
although most of the first night was spent trying differ- 
ent setup configurations, so this too produced little data. 
Typically, four sets of three spectra were obtained on each 
of the remaining nights, with exposures mostly of 1000s. 
In total 180 spectra of GP Vel were obtained. 

The echelle and CCD camera were mounted at the 
coude focus. All observations used the 31.6 lines mm -1 
echelle grating, the 158 lines mm -1 cross disperser grating, 
the 81.3 cm focal length camera and the thinned 2k x 2k 
Tektronix CCD, with a slit width of 2 arcsec. Around 70 
separate echelle orders fell onto the detector spanning a 
wavelength range 3370 - 5970 A. To shorten the readout 
time, 2x2 binning was used on the CCD for the first seven 
nights, but this was changed to 2 x 1 binning (i.e. binning 
only in the cross-dispersion direction) for the remainder 
of the run. 

The radial velocity standards HR 1829 (a G5II star) 
and HR 3694 (a K5III-IV star) were also observed on most 
nights, as well as a BOIa star (e Ori) to provide a template 
spectrum for the cross correlation. We also obtained the 
usual complement of flat fields and bias frames as well as 
thorium-argon arcs each time the telescope was moved. 

4. Data Reduction and Analysis 

The iraf software package was used to perform the re- 
duction and analysis. For each night's spectra, an aver- 
age bias frame was constructed and subtracted from the 
arc, flat and object frames. Orders were traced automati- 
cally using the GP Vel frames, as the object was relatively 
bright and its early-type spectrum does not contain many 
large absorption features which could complicate the trac- 
ing process. An optimal extraction algorithm was used to 
extract all the spectra in order to maximise the signal-to- 
noise ratio. 

A wavelength calibration solution was found for the 
whole echelle frame using the arcs. A total of 1109 
lines were identified from the thorium-argon arc spectrum 
(about 20 per order of the echelle) and the wavelength 
solution was fitted using a Chebyshev polynomial of or- 
der three in the wavelength direction plus another of or- 
der four in the spatial direction. The RMS of the fit was 
0.005A, which is less than 10% of a pixel in all orders. The 
FWHM of the individual arc lines near the centre of the 
echelle about 0.1 A. 
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Fig. 1. An average spectrum of GP Vel. Horizontal bars indicate the regions used in the cross-correlation procedure. 
Many of the sharp lines peaking downwards are weak interstellar lines; most of the upward peaking lines are remaining 
sky features. 
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At the blue end, the signal-to-noise ratio was poor 
and there was also a problem with overlapping orders. 
Consequently only 53 orders were extracted from each 
echelle spectrum, spanning the wavelength range 3820A to 
5958 A (orders 148 to 96). The spectral resolution ranged 
from 0.06 A pixel _1 at the blue end to 0.1 A pixel _1 at 
the red end. All spectra were extracted onto a common 
log-linearised scale with a resolution of 4.86 km s" 1 per 
pixel. An average spectrum is shown in Figure 1. 

The spectra of GP Vel showed essentially the same 
set of absorption lines as were seen by van Kerkwijk et 
al (1995). In particular we saw the Balmer series from 
H/3 to H9, Hei lines at 5015A, 492lA, 4713A, 447lA, 
4387A, 4143A, 412lA and 4026A, and several lines due to 
Sim (4553A, 4568A, 4575 A, 4820A, 4829A), Siiv (4089A, 
4116A), On (459lA, 4596A) and Nil (4601A). 

Before proceeding with the main cross-correlation, the 
stability of the detection system was checked by cross- 
correlating the various spectra of the radial velocity stan- 
dards HR 1829 and HR 3694 against each other. Across 
the whole run, the mean shift in these spectra was only 
0.46 km s _1 , which is negligible (< 0.1 of a pixel). As 
a further check, regions of the GP Vel spectra contain- 
ing the Call 3934A and Nal 5890/5896A interstellar lines 
were cross correlated against themselves. No trends were 
apparent in these data either. 

To extract the radial velocity curve of GP Vel, tem- 
plate spectra of e Ori were cross correlated against the 
spectra of the target. Since these stars have similar spec- 
tral types, systematic errors in this process should be min- 
imised. Regions of the spectrum containing the Hei lines 
plus the Sim and Siiv lines were used to produce the radial 
velocity curve in Figure 2 (see Appendix for the individual 
radial velocity values, relative to e Ori). The radial veloc- 
ities themselves were extracted using standard IRAF rou- 
tines by fitting a Gaussian to the cross-correlation function 
in each case. The Gaussian fits were performed using the 
full shape of the Gaussian curve down to one-quarter of 
their height. Although radial velocities were also obtained 
from the Balmer lines, these lines were in general very 
broad and so it was difficult to obtain accurate velocities. 
Furthermore, there was some evidence that the radial ve- 
locity amplitude increased with the order of the Balmer 
line. We comment on this further below, but as a result 
did not use Balmer lines in the rest of our analysis. 

5. Fitting the radial velocity curve 

Because Vela X-l is known to have a significant eccentric- 
ity, simply fitting a sinusoidal radial velocity curve to the 
cross-correlation results as a function of time would not 
be satisfactory, as an eccentric orbit will show deviations 
from a pure sinusoid. Instead it may be shown that the 
observed radial velocity is given by: 

47raisini ecosw + cos(f + u) . . 

v = 1+ p^W> 2 (10) 




True anomaly / 2ix 



Fig. 2. Radial velocity curve for Hei, Sim and Siiv lines 
showing the fitted curve (the first fit). Each cluster of ra- 
dial velocity values represents a single night's data. The 
fit has a reduced chi-squared of 3.84. 

where a\ is the semi-major axis of the orbit of GP Vel, 
e is the orbital eccentricity, P is the orbital period, v is 
the true anomaly (i.e. the angle between the major axis 
and a line from the star to the focus of the ellipse), u is 
the periastron angle (i.e. the angle between a line from 
the centre of the orbit to periastron and a line from the 
centre of the orbit to the ascending node), i is the angle 
between the normal to the plane of the orbit and the line- 
of-sight and 7 is the radial velocity of the centre of mass 
of the binary system. The true anomaly is related to the 
eccentric anomaly, E, by: 

'*■■(!) 

where E is the angle between the major axis of the 
ellipse and the line joining the position of the object and 
the centre of the ellipse. In turn, E can be related to M, 
the mean anomaly, by Kepler's Equation: 

2tt 

E-esinE = M= — (t - T ) (12) 

where the right hand side of this equation is simply 
the orbital phase with To the time of periastron passage. 
This equation can be solved numerically, or using Bessel 
functions. 

The value of T was obtained from the value of T v / 2 > 
the epoch of 90° mean orbital longitude, derived by 
Bildstcn et al (1997) from BATSE pulse timing data. They 
give T„ /2 = MJD 48895.2186 ± 0.0012, from which 

To = T w/2 + V 2?r 21 (13) 

where w — 152.59°±0.92°. Hence we calculate the time 
of periastron passage as T = MJD 48896.777 ± 0.009. 



6 



Quaintrcll, Norton, Ash et al: Neutron star mass in Vela X-l 




115 120 125 130 135 



HJD - 2450000 

Fig. 3. Residuals to the radial velocity curve fit in Figure 
2, plotted against time. (The first fit.) Overlaid is a best-fit 
sinusoid with a period of 2.18 d. 

5.1. The first fit 

Times, t, were assigned to each radial velocity measure- 
ment and Equation 12 was solved using a numerical grid 
to calculate the corresponding eccentric anomaly. The true 
anomaly was then found using Equation 11, and Figure 2 
shows a plot of radial velocity against true anomaly, with 
the best fit curve according to Equation 10. The reduced 
chi-squared of the fit is x? = 3.84. Scaling the error bars 
by a factor of 1.96 to reduce x? to unity, gives the am- 
plitude of the fitted curve as K = 21.4 ± 0.5 km s _1 . 
However, we note that, since the use of chi-squared as- 
sumes that the errors on all the points are uncorrelated, 
the uncertainty here is likely to be a gross under-estimate. 
This will hereafter be referred to as the 'first fit'. 

Figure 3 shows the residuals to the radial velocity curve 
in the first fit, plotted against time. It is clear that there 
are trends apparent in these residuals from night to night 
and a Fourier analysis shows that the dominant signals 
are at periods of 9 ± 1 d and 2.18 ± 0.04 d (Figure 4). The 
9 d signal present in the power spectrum reflects the fact 
that fixing the zero phase of the radial velocity curve, as 
implied by the X-ray data of Bildsten et al (1997) does 
not provide the best fit to the data. We suggest this effect 
may be responsible for the 'phase-locked' deviations re- 
vealed by Barziv et al (2001) too. The 2.18 d modulation 
appears to be relatively stable throughout our two orbits 
of observations, with an amplitude of around 5 km s _1 , as 
shown in Figure 3. 

5.2. The second fit 

In order to take account of the correlated residuals from 
the first fit, another fit to the radial velocity data was per- 
formed this time with four free parameters: the 7 velocity 
and the K Q value as before, plus the amplitude and phase 




Fig. 4. Power spectra of the residuals to the radial veloc- 
ity curve fits. The solid curve is the power spectrum of 
the data in Figure 3 (the residuals after the first fit) and 
the dotted curve (offset in the — ve direction by one unit 
for clarity) is the power spectrum of the data in Figure 
6 (the residuals after the second fit). The highest peaks 
in the power spectrum of the residuals after the first fit 
correspond to frequencies of 0.111 d _1 and 0.459 d _1 or 
periods of 9.0 d and 2.18 d. In each case, window function 
peaks due to the sampling of the lightcurve have been re- 
moved using a 1-dimensional clean (program courtesy H. 
Lehto) . 

of a sinusoid with a period of 2.18 d. The K part of the 
fit is sinusoidal as a function of true anomaly, whilst the 
2.18 d part of the fit is sinusoidal purely as a function of 
time. The radial velocity fit according to this prescription, 
referred to as the second fit, is shown in Figure 5 and has 
a reduced chi-squared of 2.40. The improvement over the 
first fit is therefore clearly apparent. Scaling the error bars 
by a factor of 1.55 to reduce x? to unity, gives the ampli- 
tude of the fitted curve as K Q = 22.4 ±0.5 km s _1 and the 
amplitude of the 2.18 d modulation as 5.4±0.5 km s _1 . As 
with the first fit, we note that the uncertainties here are 
also likely to be under-estimates. The remaining residuals 
on the second fit are shown in Figure 6, and their power 
spectrum is indicated by the dotted line in Figure 4. It can 
be seen that some of the correlated structure in the resid- 
uals has been removed by this procedure, although the 
peak corresponding to a period of ~ 9 d is still present 
indicating once again that the zero phase according to the 
X-ray observations does not provide the best description. 

5.3. Fits to phase binned means 

To test whether the deviations from the radial velocity 
curve displayed by each group of data points are corre- 
lated, we subtracted the effect of the 2.18 d period ac- 
cording to the parameters found for it in the second fit, 
and rebinned the resulting radial velocities into nine phase 
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Fig. 5. Radial velocity curve for Hei, Sim and Siiv lines 
showing the fitted curve (the second fit) after removal of 
the 2.18 hour signal. Each cluster of radial velocity values 
represents a single night's data. The fit has a reduced chi- 
squared of 2.40. 
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Fig. 6. Residuals to the radial velocity curve fit in Figure 
5, plotted against time. (The second fit.) 

bins (each ~ 1 day) . The radial velocity fit to these phase- 
binned means is shown in Figure 7 and has a reduced chi- 
squared of 28.3. The errors on the mean data points are 
calculated as the standard error on the mean from the in- 
dividual data points that are averaged in each case. The 
fact that the reduced chi-squared is significantly greater 
than unity suggests that the deviations are indeed cor- 
related. Scaling the error bars here by a factor of 5.3 to 
reduce x? to unity, gives the amplitude of the fitted curve 
as K = 22.6 ± 1.5 km s" 1 . 

Figure 7 clearly shows that a better fit would be ob- 
tained if the zero phase of the radial velocity is allowed 
as a free parameter. This reflects the fact that a ~ 9 d 
period is present in the power spectra of the residuals af- 




0.5 1 



True anomaly / 27T 

Fig. 7. Radial velocity curve for Hei, Sim and Siiv lines 
binned into nine phase bins. Each radial velocity value 
represents two or more night's data. The fit to the phase 
binned means has a reduced chi-squared of 28.3 and is 
shown by the solid curve. A second fit to the phase binned 
means, allowing the zero phase of the radial velocity as a 
free parameter, is shown by the dashed line and has a 
reduced chi-squared of 6.0. 

ter both the first and second fits (Figure 4). Allowing the 
zero phase to vary results in the second fit to the phase 
binned data shown by the dashed line in Figure 7. This 
time the reduced chi-squared is 6.0. Scaling the error bars 
here by a factor of 2.44 to reduce Xr to unity, gives the 
amplitude of the fitted curve as 21.2 ± 0.7 km s _1 . The 
shift of the minimum radial velocity with respect to that 
in the first fit to the pase binned means is a true anomaly 
of Av = 0.033±0.007, corresponding to about 7 hours. We 
note that this shift in zero phase is significantly greater 
than the uncertainty implied by extrapolating the X-ray 
ephemeris of Bildsten et al (1997) to the epoch of our ob- 
servations which is only ~ 0.0006 of a cycle. We conclude 
that the orbital phase locked deviations remaining in the 
radial velocities after the first fit to means can therefore 
mostly be accounted for by a simple phase shift in the 
radial velocity curve. 

6. Calculation of system parameters 

The masses of the two stellar components were calculated 
according to the method described in Section 1. In order 
to propagate the uncertainties, a Monte Carlo approach 
was adopted, calculating 10 6 solutions with input values 
drawn from a Gaussian distribution within the respective 
la uncertainties of each parameter. The final values were 
then calculated as the mean result of these 10 6 solutions, 
with the lcr uncertainties given by the rms deviation in 
each case. 

The procedure was carried out using K Q = 22.6 ± 
1.5 km s _1 from the first fit to means and the results 
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of these calculations are shown in Table 2. We used the 
value of K x = 278.1 ± 0.3 km s" 1 (Bildsten et al 1997), to 
calculate q using Equation 1 and a value of 33°. 8 ± 1.3° 
was adopted for the eclipse half-angle, 9 e from Watson & 
Griffiths (1977). We note that this value was also used by 
Wilson & Terrel (1998) for their unified analysis of Vela X- 

1 and appears to be the most reliable. A value of 0.67±0.04 
was adopted for the co-rotation factor, 0, obtained from 
comparison of GP Vel spectra with model line profiles by 
Zuiderwijk (1995). We note that Barziv et al (2001) also 
used the wsini value derived by Zuiderwijk (1985) and 
obtained a similar value, namely il(= f co ) = 0.69 ± 0.08. 

Values for (3 and i cannot be obtained from indepen- 
dent observations. However, since the supergiant is un- 
likely to be overfilling its Roche lobe, an upper limit to 
the Roche lobe filling factor is (3 < 1.0; and an upper limit 
to the inclination is clearly i < 90°. A lower limit to (3 is 
found by setting i — 90° , and a lower limit to i is found by 
setting (3 — 1. These limits are indicated in Table 2. We 
note though that (3 is unlikely to have a fixed value, as the 
size and shape of the Roche lobes of both components will 
vary as the separation between the stars varies as they or- 
bit each other eccentrically. Indeed, it has been suggested 
that GP Vel is likely to fill its Roche lobe at periastron 
(Zuiderwijk 1995). However the mass ratio implies that 
Roche lobe overflow would be dynamically unstable and 
lead to a common-envelope phase, so the limit of (3 = 1 
is probably implausible. We also note that Barziv et al 
(2001) took a somewhat different approach and assumed 
a filling factor at periastron in the range (3 = 0.9 — 1.0, 
thus leading to a minimum value for the inclination of 73° , 
consistent with our analysis. 

Given these two limits on each of (3 and i, Equations 

2 and 3 then allow us to calculate the masses of the two 
components for combinations of i and (3 between the two 
extremes. The limiting values for the masses of the two 
stars are shown in Table 2. Values between these extremes, 
as a function of inclination angle, are shown in Figure 8. 
With these two limiting situations, we also calculate the 
semi major axis a of the orbit from Kepler's law, and this 
is then used to calculate the separation of the stars at peri- 
astron according to a' = a(l — e). The Roche lobe radius at 
periastron Rl is then determined from the value of Rl/ci' 
at each extreme, and the radius of the companion star R a 
is then calculated as j3 x i?L in each case. These values 
too are shown in Table 2. For each parameter, the actual 
value may be anywhere within the range encompassed by 
the two limits. 

7. Tidally induced non-radial oscillations 

The strongest peaks in the power spectrum of the residuals 
after the first fit, at 0.111 d _1 and 0.457 d _1 , are sugges- 
tively close to the orbital frequency and 4x the orbital 
frequency (f or b = 0.11155 d -1 ) respectively. This leads us 
to speculate that they may both represent tidally excited 
non-radial oscillation modes (see also Willems & Aerts 
2002, Handler et al. 2002). 



Table 2. System parameters for Vela X-l / GP Vel. 
The value for K a is that resulting from fitting the 
phase bin means without a phase shift. 
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<i 


0.081 ± 0.005 
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1.000 


0.89 ±0.03 




i / deg 


70.1 ±2.6 


90.0 




M x / M Q 


2.27 ±0.17 


1.88 ±0.13 




M / M Q 


27.9 ±1.3 


23.1 ±0.2 




a' at periastron / Rq 


51.4 ±0.8 


48.3 ±0.3 




i?L at periastron / 


32.1 ±0.6 


30.2 ±0.2 




Ro 1 R© 


32.1 ±0.6 


26.8 ±0.9 





[1] Bildsten et al 1997 
[2] Watson & Griffiths 1977 
[3] Zuiderwijk 1995 
[4] this paper 



From a theoretical point of view, the tidal action ex- 
erted by one binary component on the other is governed 
by the tide-generating potential which can be expanded in 
a Fourier series in terms of multiples of the mean motion 
n = 2 7r/ or b (e.g. Polfliet & Smeyers 1990). The oscilla- 
tions that are most likely to be excited are those associ- 
ated with the lower-order harmonics of the Fourier series. 
In the particular case of GP Vel, the dominant contribu- 
tions to the Fourier series of the tide-generating potential 
are associated with the first ten harmonics. This supports 
the possibility that the oscillation associated with the fre- 
quencies « /orb and « 4 / or b may be induced by the tidal 
action of the neutron star. 

An in-depth analysis of the possibility that a tidally 
excited oscillation mode exists in GP Vel requires an ac- 
curate knowledge of the stellar and orbital parameters. 
In particular, the internal structure of the B-type com- 
panion to the neutron star must be known to derive the 
spectrum of the eigenfrequencies present in the star. The 
radii quoted in Table 2 indicate that the B-type star might 
be in the Hertzsprung-gap. An additional complication 
arises from the unknown evolutionary history of the bi- 
nary as the B-type star may be 'polluted' by mass-transfer 
episodes prior to the formation of the neutron star. In view 
of these complications, a detailed study of the oscillation 
spectrum of GP Vel and its possible use in an asteroseis- 
mological study is beyond the scope of our present inves- 
tigation. 
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Fig. 8. The masses of the neutron star (upper panel) and 
companion star (lower panel) in solar units, as a function 
of inclination angle. The lower limit to the inclination an- 
gle is set by the condition that the companion star cannot 
over-fill its Roche lobe. Error bars on the masses at the 
extreme limits are determined as explained in the text. 

8. Discussion 

Our final value for K of 22.6±1.5 km s" 1 is slightly higher 
than, but comparable to, those quoted by other recent 
authors, noted in Section 2. Our result for the neutron 
star mass is compatible with van Kerkwijk et al's (1995) 
mass value of M x = I.9+O5M0, and that of Barziv et al 
(2001), namely M x = 1.86±O.32M . We therefore support 
the view that the neutron star in Vela X-l is more massive 
than the canonical value of 1.4M0. However, there are a 
number of potential sources of systematic error, as noted 
in our earlier paper concerning the mass of the neutron 
star in Cen X-3 (Ash et al 1999), some of which will be 
examined below. 

8.1. Stellar wind 

As GP Vel is an early-type star, it has a significant stel- 
lar wind, and a strong and variable stellar wind can have 



disastrous consequences when making radial velocity mea- 
surements, especially over an extended period of time. 
Barziv et al (2001) made a detailed study of the effects 
of winds on their observations and showed that deviations 
in the Balmer lines in particular indicated the presence 
of a photo-ionization wake. We note that the supergiant 
GP Vel is of later type (B0) than that of Cen X-3 (06-7), 
and thus has a weaker wind. Additionally, as our observa- 
tions were made over a relatively short time period, the 
potential effect of long-term wind variation will be less 
than if the observations were spread over a number of 
years, as in the case of our Cen X-3 campaign, or over a 
year, as in the case of Barziv et al's Vela X-l campaign. In 
other words, although the effect of the supergiant 's stel- 
lar wind cannot be discounted, its effect should be less 
than its effect on our Cen X-3 results and less than on the 
results of Barziv et al (2001). 

The apparent increase in amplitude of the radial veloc- 
ity curve seen with increasing order in the Balmer series 
lines mentioned earlier, may also be an indication of a 
stellar wind. Conceivably, the different order lines have 
different contributions from the wind material as it moves 
further out, and this could produce the correlation seen. 
Since we have not used the Balmer lines in our analysis 
however, this effect does not influence our final result. 

8.2. X-ray heating 

The effect of X-ray heating of the surface of GP Vel may be 
compared with the corresponding effect in Cen X-3. The 
ratio of intrinsic to irradiated fluxes may be estimated for 
each system as (L / AttR 2 ) / (L^/ Air(a - R ) 2 ). The RXTE 
ASM lightcurves show that both sources have similar X- 
ray fluxes, but Cen X-3 is at over three times the distance 
of Vela X-l (Sadakane et al 1985; Humphreys & Whelan 
1975) so its X-ray luminosity is at least a factor of ten 
higher. The luminosities of the supergiant stars in both 
systems are similar, but the radius of the star in Cen X-3 
and the separation of the supergiant and neutron star are 
both around three times smaller in Cen X-3 (Ash et al 
1999) than in Vela X-l. Hence the effect of X-ray heating 
in Vela X-l is of order ten times less than in Cen X-3. 
We were able to rule out significant X-ray heating in the 
case of Cen X-3 (Ash et al 1999), and the effect should 
therefore be even less in this case. 

8.3. Deviations from Keplerian radial velocity curve 

Clearly, another source of uncertainty in the case of 
GP Vel arises from the observed deviations from a pure 
Keplerian radial velocity curve, van Kerkwijk et al (1995) 
suggested that these deviations do not last longer than a 
single night, although our results (Figure 3) suggest that 
they may last for up to two system orbits. They suggested 
that the deviations are the result of short-lived, high-order 
pulsations of the photosphere and we confirm that the fre- 
quencies appear to be harmonically related to the orbital 
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frequency. It is possible therefore that they do represent 
tidally induced oscillations. As well as introducing another 
source of error into our attempts to measure the radial 
velocity of the centre of mass of the supergiant, the "wob- 
bles" cause another problem. If the shape of the star is 
constantly changing, this will affect the extent to which 
the star fills its Roche lobe. Additionally, the size and 
shape of the Roche lobe itself will vary with phase, as the 
two components orbit each other in eccentric orbits. We 
believe that we have partially accounted for this effect by 
removing the 2.18 d modulation shown in Figure 3. 

However, as Figure 7 demonstrates, there is clearly a 
remaining effect which may be characterised by an addi- 
tional radial velocity variation at the orbital period. This 
manifests itself as a phase shift in the radial velocity curve 
with respect to the X-ray orbital ephemeris. The measured 
amplitude of the radial velocity curve is therefore a com- 
bination of the orbital motion of GP Vel and an additional 
component that may be due to a non-radial pulsation in- 
duced by the eccentric orbit of the neutron star. There 
may even be another component of this non-radial pulsa- 
tion at the orbital period that is in phase with the orbital 
motion and so undetectable. Our K a value, and presum- 
ably all previous measurements too, are therefore subject 
to an unknown systematic error, and the stellar masses 
calculated from them are also uncertain because of this. 

Another possibility for the origin of the extra velocity 
component at the orbital period is large-scale motion re- 
quired to keep the tidal bulge pointing in the direction of 
the neutron star given the asynchronous stellar rotation 
(fi = 0.67). However, such a tidal bulge would extend in 
both directions - towards and away from the neutron star 
- and so the effect may cancel out. 

As a further result, our identification of these depar- 
tures from the radial velocity curve in GP Vel suggest that 
it may be the third recent candidate to exhibit tidally 
induced non-radial oscillations, after HD 177863 (De Cat 
2001, Willems & Aerts 2002), and HD 209295 (Handler et 
al. 2002). 

8.4. Neutron star mass 

The minimum neutron star mass consistent with our re- 
sults is around 1.74 M Q . A recent review of neutron 
star equations of state by Lattimer & Prakash (2000) 
shows that such a mass would exclude the schematic 
potential models of Prakash, Ainsworth and Latimer, 
and the various models based on a field theoretical ap- 
proach by Glendenning & Moszkowski, by Glendenning 
and Schaffner-Bielich, and by Prakash, Cooke & Lattimer. 
Hence the 'softer' equations of state that incorporate com- 
ponents such as hyperons, kaons and quark matter are 
ruled out. 

In order to be consistent with a neutron star mass of 
~ 1.4 Mq, the radial velocity amplitude of GP Vel would 
have to be K D ~ 18 km s -1 , implying that the additional 



velocity component due to the non-radial pulsation has an 
amplitude of around ~ 4.5 km s _1 . 

We conclude by noting that there are two large sources 
of uncertainty in deriving the masses of the stellar compo- 
nents. The first is that of the value j3 (the ratio of the ra- 
dius of the companion star to its Roche lobe radius). This 
in turn leads to a large uncertainty in sin i (and hence the 
masses) and is unlikely to be reduced even with optical 
spectroscopic observations of greater signal-to-noise. The 
second source of uncertainty is the additional radial veloc- 
ity component at the orbital period that we have discov- 
ered and represented by a phase shift. Removing the effect 
of such a velocity component is likely to require detailed 
modelling of tidally excited oscillation modes in GP Vel 
which in turn requires an accurate knowledge of the stel- 
lar parameters. Consequently, the mass measurement pre- 
sented here may represent the best that is achievable. 
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Table 3. The radial velocities of GP Vel relative to that of e Ori 



HJD - 2450000 


RV / km s" 


i 


HJD - 2450000 


RV / km s" 


L 


HJD - 2450000 


RV / km s" 


L 


115.9492 


-39.978 ± 2 
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124.1798 


-52.602 ±2 


338 


129.9264 


-70.958 ± 2 


530 


116.1643 


-33.359 ± 2 


026 


124.1944 


-50.835 ± 2 


627 


129.9393 


-70.168 ±2 


452 


116.1768 


-32.720 ± 2 


218 


124.2083 


-51.151 ±3 


466 


129.9995 


-70.815 ±2 


380 


116.1905 


-32.454 ± 2 


090 


124.9819 


-30.123 ±2 


138 


130.0124 


-71.690 ±2 


064 


116.9187 


-31.398 ±2 


249 


124.9944 


-29.702 ± 2 


094 


130.0251 


-70.973 ± 2 


494 


116.9314 


-30.109 ±2 


933 


125.0078 


-29.527 ±2 


370 


130.0515 


-71.616 ± 2 


018 


116.9468 


-30.908 ± 3 


410 


125.1431 


-27.910 ±2 


830 


130.0650 


-71.685 ±2 


830 


117.0059 


-34.455 ± 5 


798 


125.1561 


-28.457 ± 3 


213 


130.0780 


-72.801 ±3 


691 


117.9872 


-46.720 ± 3 


585 


125.1691 


-26.982 ± 2 


189 


130.1179 


-70.146 ± 2 


440 


118.0148 


-47.951 ± 6 


251 


125.9233 


-34.387 ± 2 


657 


130.1314 


-69.869 ± 2 


277 


118.0710 


-48.042 ± 3 


267 


125.9373 


-34.510 ±2 


742 


130.1447 


-69.998 ± 2 


608 


118.0877 


-51.834 ±5 


119 


125.9506 


-34.692 ± 2 


958 


130.1865 


-69.660 ± 3 


002 


118.1010 


-53.321 ± 7 


217 


126.0339 


-33.610 ±2 


408 


130.1996 


-70.063 ± 2 


521 


118.1454 


-50.003 ± 2 


466 


126.0469 


-33.979 ± 2 


359 


130.2126 


-70.230 ±3 


305 


118.1585 


-48.447 ± 5 
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-56.218 ±2 


664 


120.1824 


-75.016 ±2 


495 


128.0551 


-62.303 ±4 


473 


133.0305 


-54.996 ± 2 


741 


120.1962 


-74.390 ± 2 


519 


128.0698 


-61.858 ±4 


231 


133.0435 


-53.878 ± 2 


518 


120.2103 


-74.013 ± 2 


769 


128.0829 


-61.410 ±4 


222 


133.0570 


-54.587 ± 2 


539 


120.2417 


-72.921 ±2 


729 


128.1450 


-64.184 ±4 


351 


133.9029 


-35.155 ±3 


759 


120.2543 


-72.066 ± 2 


847 


128.1585 


-61.980 ±3 


912 


133.9158 


-36.872 ± 3 


280 


120.2669 


-74.260 ± 3 


256 


128.1715 


-60.965 ± 4 


288 


133.9288 


-37.601 ±3 


099 


120.9177 


-69.011 ±2 


892 


128.1846 


-65.137 ±3 


990 


133.9773 


-36.505 ± 3 


376 


120.9298 


-67.862 ± 2 


769 


128.2014 


-63.062 ± 3 


888 


133.9903 


-35.226 ±3 


121 


120.9430 


-68.414 ±2 


985 


128.2145 


-63.871 ±4 


560 


134.0064 


-36.135 ±3 


172 


120.9659 


-67.933 ± 2 


808 


128.2275 


-63.709 ± 4 


447 


134.0861 


-37.570 ± 3 


744 


120.9789 


-68.716 ±2 


809 


128.9267 


-69.557 ± 3 


368 


134.9085 


-28.154 ±3 


774 


120.9929 


-67.612 ±2 


848 


128.9396 


-69.768 ± 3 


061 


134.9215 


-29.949 ± 2 


582 


121.0343 


-67.387 ± 2 


906 


128.9526 


-70.041 ± 3 


134 


134.9346 


-29.342 ± 2 


227 


122.1584 


-68.859 ± 5 


309 


129.0007 


-68.518 ±3 


256 


134.9852 


-28.680 ± 2 


028 


122.1803 


-76.832 ± 5 


303 


129.0142 


-70.118 ±3 


534 


134.9982 


-29.925 ± 2 


816 


122.1934 


-70.372 ± 3 


793 


129.0272 


-68.765 ± 3 


262 


135.0146 


-29.764 ± 2 


670 


123.9486 


-55.418 ±2 


464 


129.1112 


-68.900 ± 2 


794 


135.0596 


-30.751 ±3 


978 


123.9680 


-53.818 ±2 


188 


129.1244 


-68.574 ± 2 


835 









